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Defects often noticeably influence the electrical and op- 
tical properties of a material by introducing defect states 
into the band gap. Reaching a microscopic understand- 
ing of the physical and chemical properties of defects in 
solids has long been a goal of first-principles electronic 
structure methods. Probably the most widespread theo- 
retical method in this realm today is density functional 
theory (DFT) in the local-density (LDA) and general- 
ized gradient approximation (GGA), but certain intrin- 
sic deficiencies limit their predictive power. Artificial 
self-interaction and the absence of the derivative discon- 
tinuity in the exchange-correlation potential [1J present 
the most notable deficiencies in this context. They give, 
amongst other things, rise to the band-gap problem - the 
fact that the band gap in LDA and GGA underestimates 
the quasiparticle gap [HE]- In this Letter we show that 
the band-gap problem in LDA / GGA not only affects the 
reliable computation of defect levels, but in certain cases 
(e.g. filled defect states in the band gap) also that of for- 
mation energies. We present a formalism for calculating 
formation energies of defects in solids that combines LDA 
with quasiparticle energy calculations in the GqWq ap- 
proximation [3j to reduce the self-interaction error and to 
overcome the band-gap problem. In some cases a heuris- 
tic "scissor operator" approach may approximately cor- 
rect the problem. However, particularly when the exper- 
imental answer is unknown, a more accurate method is 
needed. 

We illustrate our approach with the example of a self- 
intersitital in silicon (Si;), a defect of high technological 
relevance |H |5J |B] . In the neutral charge state the Si; has 
several stable and metastable atomic configurations [Tj |S] 
(see Fig. [I}, in all of which two electrons occupy a de- 
fect level in the band gap. The LDA formation energies 
of all these configurations are underestimated severely 
(by ~1.5 eV) compared to diffusion Monte Carlo (DMC) 
calculations [5J [TO]. However, no insight into this dis- 




FIG. 1: (Color online) a) Split<110>, b) hexagonal, c) C 3v 
and d) tetrahedral configuration of the Sij . Defect atoms are 
shown in red and nearest neighbours in grey. The middle 
panel depicts the formation of the neutral Sii from the 2+ 
charge state. A+ and A2+ are short for the electron affinities 
j4(+,Ro) and A(2+, R2+) (see text), respectively, and R 9 
denotes the atomic positions in charge state q. 



crepancy is provided by the DMC calculations. 

In our formalism the formation of the neutral defect is 
expressed as successive charging of its 2+ charge state, 
for which the defect level is unoccupied. This allows us 
to decompose the formation energy into that of the 2+ 
state (E' (2+)), a lattice and an electron addition part. 
This decomposition is not only insightful for analyzing 
the underestimation of the LDA formation energy, but 
also permits us to apply the most suitable method for 
each type of contribution [11] . For the lattice part we re- 
tain the LDA and argue that the relaxation energies and 
E' (2+) are not as strongly affected by the deficiencies of 
the LDA as in the positive and the neutral case since the 
defect level in the band gap is unoccupied. For the elec- 
tron affinities, on the other hand, we employ Hedin's GW 
method |3] . Since self-consistency in GW is still discussed 
controversially [2 we obtain the quasiparticle corrections 
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to the LDA Kohn-Sham energies from first order pertur- 
bation theory (Go Wo), which is currently the method 
of choice for computing quasiparticle band structures in 
solids [2J [TJ] . While not completely self-interaction free 
P2] Go Wo significantly reduces the self-interaction error. 
With this combined approach the formation energy in the 
neutral charge state increases by ~1.1 eV compared to 
the LDA. Recent DFT calculations employing the HSE 
hybrid functional, which also significantly reduces the 
self-interaction error, yield a similar improvement |10j 
and lend further substance to this notion. Moreover, the 
GoWo-corrected charge transition levels agree well with 
recent experimental measurements [4]. 

We will present our combined DFT+GoWo approach 
for the example of the Si;, but it can easily be generalized 
to defects in other materials. An additional silicon atom 
in an interstitial site can adopt different configurations 
with similar formation energies (cf. Fig. [TJ. In the tetra- 
hedral (tet) geometry the extra silicon atom gives rise 
to an ai and a threefold degenerate state. The lat- 
ter is empty and in resonance with the conduction band. 
The partial occupation of the degenerate ti state triggers 
a Jahn- Teller distortion along the <lll>-axis into a ge- 
ometry with Cs v symmetry, also referred to as "displaced 
hexagonal structure" in previous studies |14j . The addi- 
tion of a second electron displaces the atom further in 
the <lll>-direction stabilizing the neutral charge state. 
This sequence is illustrated in Fig. [2] Moving the inter- 
stitial atom along the <lll>-direction into the center 
of a six-membered ring (hexagonal (hex) configuration) 
lowers the neutral state further in energy. It reaches its 
lowest position in the split <110> configuration, where 
the added atom and a host atom share an interstitial site 
oriented in the <110>-direction. 

For the 2+ charge state the tetrahedral is the most 
stable configuration [5] (see also Table [I]) and we will use 
it as a starting point for building our scheme. The pos- 
itive charge state is then formed by adding one electron 
as depicted in steps 1 and 2 in Fig. [l] Mathematically 
this can be expressed by starting from the expression for 
the formation energy in the positive charge state 

E f D (+,eF) = E(+,TL°)-E ref + e F . (1) 

E(q, Ry ) is the total energy in charge state q and atomic 
positions Ry of defect configuration D in charge state 
q' . E re f is a suitably chosen reference system, here bulk 
silicon, and tp the Fermi level of the electrons referenced 
to the top of the valence band. Adding and substracting 
first £:(+,RJf_*) and then £'(2+,R|!f) leads to 

E f D (+, e P ) =A(+, R?,B£f) + A(2+, R£*) 

+ E{ et (2+,e F = 0) + e F . (2) 

The energy difference £7(+,R|^) - E(2+,R t 2 °l) defines 
the vertical electron affinity A(2+, R2+) of the 2+ state 



(in its tetrahedral configuration), step 1 in Fig. [T| 
referenced to the top of the valence band, whereas 
i?(+,R+) — £ , (+,R|+) gives the subsequent relaxation 
energy A(+, R^, R2+) in the positive charge state (step 
2). _ 

Similarly the neutral charge state emerges from the 
positive one by addition of an electron. Mathematically 
we again achieve this by adding and substracting first 
E(+, Rq) and then E(+, R+ ) to and from the expression 
for the neutral formation energy Ejj(0, ep) = E(0, Hq) — 

E re f'- 

E f D (0, e F ) =A(+, Rf ) + A(+, R?,R?) 

+ Ef(+,e F =0) . (3) 

Again A(+,Rq) denotes a vertical electron affinity 
E(0,Rg) - E(+,Rq) (step 4) and A(+, R^, R^) = 
E(+,RV) - E(+,TL%) the relaxation energy from the 
neutral to the positive geometry in the positive charge 
state (step 3) . An expression for the negative charge state 
can be obtained completely analogously once Ejj(0, ep = 
0) has been computed. 

The decomposition in Eq. Q and ([3| is not only ap- 
pealing from an intuitive point of view, but also groups 
the required total-energy differences into two categories: 
lattice contributions in a fixed charge state and electron 
addition energies at fixed geometry. This permits us to 
go beyond a pure DFT description in an easy fashion 
by employing the most suitable method for each type of 
contribution Since we expect relaxation energies in 
the same charge state to be given reliably by LDA we 
retain DFT for the lattice part. For the electron ad- 
dition energies, i.e., changes in charge state, which are 
typically problematic in LDA, we instead resort to the 
G Wo method. 

The last remaining quantity to be assigned is 
E{ et (2+,ep = 0), which we compute in the LDA. Un- 
like for the neutral state, the absence of DMC reference 
data unfortunately does not permit an assessment of the 
LDA error in this case. However, since the conduction- 
band-derived defect levels are unoccupied the effect of 
the self-interaction and the band-gap error on the forma- 
tion energy should be small. We therefore expect LDA 
to be more reliable for the tetrahedral 2+ state than for 
the neutral or the positive states. 

The LDA calculations in the present work have been 
performed with the plane-wave, pseudopotential code 
S/PHI/nX 15]. 64-atom supercells were used through- 
out, unless otherwise noted. To remove the contributions 
arising from the homogeneous compensation charge den- 
sity that is added to charged supercell calculations we 
have performed calculations for supercells with 64, 216 
and 512 atoms. In these the interstitial atom was placed 
in the tetrahedral (2+) and the C^ v (+) position of a 
perfect (undistored) Si lattice. Fitting the formation en- 
ergies up to cubic order in the inverse cell length and 
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TABLE I: GoWo vertical electron affinities for different Sii configurations D and LDA relaxation energies. A( — 1 R^^Rq > ^) is 
-0.028 eV for the split<110> and A(2+, R2+) amounts to 1.26 eV in Go Wo. The tetrahedral configuration is taken as the 2+ 
state of the C3i, . Corrections for charged supercells (see text) have been added. All values are given in eV. 
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FIG. 2: (Color online) Vertical electron affinities for different 
configurations of the Sii : LDA Slater transition states (blue) 
and GqWo quasiparticle energies (red). 



extrapolating to infinite length we obtain corrections to 
the 64-atom cell of 0.17 eV and 0.04 eV for the 2+ and 
+ state, respectively. Our extrapolated formation energy 
for the unrelaxed tetrahedral 2+ state of 3.19 eV agrees 
well the 3.31 eV obtained by Wright and Modinc for a 
slightly larger lattice constant [16 . With this correc- 
tion the formation energy of the relaxed tetrahedral 2+ 
configuration amounts to 2.65 eV. Applying a recently 
developed improved correction scheme [17] yields a cor- 
rected value of 2.66 eV, in excellent agreement with our 
extrapolated value. 

For the Go Wo calculations [18] we have employed the 
GqWq space-time code gwst []30 [201 121] • For compu- 
tational convenience we calculate the electron affinity of 
positive charge states (A(+, Rq)) by their inverse pro- 
cess, the electron removal from the neutral state, since 
no spin polarization or partially filled defect states are 
encountered then. Separate GqWq calculations for bulk 
silicon yield a band gap of 1.27 eV in good agreement 
with previous pseudopotential Go Wo calculations |12j . 

The computed vertical electron affinities are shown in 
Figure [2] For comparison the LDA affinities calculated 
as Slater transition states [55] at half occupation have 
been included. The Go Wo corrections for the +/0 state 
are similar for the three configurations and relatively 
small (~0.2 eV). For states that in the LDA are in reso- 
nance with the conduction band, however, the Go Wo cor- 
rections are much more pronounced. Since these states 
have a contribution from delocalized conduction-band 



states the derealization error of the LDA [23] leads to 
a breakdown of Slater transition state theory. The re- 
sulting severe underestimation of the vertical affinities 
is akin to the band-gap problem. In LDA the band 
gap Eg = I — A, where I is the ionization potential 
and A the electron affinity, is underestimated regard- 
less of whether / and A are calculated as total energy 
differences or by Kohn-Sham eigenvalues, because the 
exchange-correlation functional is a continuous function 
of the electron density and therefore does not exhibit a 
derivative discontinuity. Many-body perturbation theoy 
in the GW approach, on the other hand, does not suffer 
from this problem. 

Having identified the relevant electron affinities we can 
now return to the formation energies in Eqs. Q and ((3}. 
Table [TJshows that already upon adding the first electron 
to the 2+ state we observe a large correction (~0.9 eV) 
for the formation energy of the positive state. This er- 
ror subsequently carries over to the neutral charge state, 
and adding the second electron incurs a further increase. 
The GoWrCorrected formation energies are now on aver- 
age 1.1 eV larger than in the LDA. Since the quasiparticle 
shift of the empty defect state in the split<110> config- 
uration is smaller than the band-gap opening the state 
is moved into the band gap (-4(0, Rg plit )=l.l eV). As a 
result the negative charge state becomes stable in Go Wo, 
which is not the case in LDA, and has a formation energy 
of 5.53 eV. 

For the neutral charge state our Go Wo corrected for- 
mation energies compare well with recent DMC calcula- 
tions that find an average increase of ~1.5 eV (with a 
statitistical error bar of ±0.09 eV) and DFT HSE hybrid 
functional calculations that significantly reduce the self- 
interaction error and yield an average increase of ~ 1.2 eV 
[10 . Earlier DMC calculations give a larger average in- 
crease of 1.7 eV compared to the LDA, but also a much 
larger statistical error bar (±0.48 eV) [S]. Assuming a 
migration barrier of ~0.2 eV [9] our computed activa- 
tion enthalpy (formation energy + migration barrier) of 
~4.7 eV for the neutral split <110> interstitial is also in 
very good agreement with the experimentally determined 
value of 4.95 eV [H]. 

Finally we address the stability of the different defect 
configurations when the Fermi energy is varied through- 
out the band gap (cf Fig. [3]). For clarity this is shown 
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FIG. 3: (Color online) Formation energies (E^) as a function 
of Fermi energy in LDA (left) and Go Wo (right). The lower 
panels show the split<110> as representative configuration 
and the upper the configuration with the lowest energy for a 
give Fermi level. The dotted line marks the LDA band gap. 



only for the split<110> configuration (lower panels) and 
the configuration with lowest energy at a given Fermi 
level (upper panels). The situation for the hex and 
C% v configurations is qualitatively similar to that of the 
split<110>. If each configuration is considered sepa- 
rately, the formation energy diagram looks startlingly 
different in LDA and Go Wo. Since LDA underestimates 
the formation energies of the + and state relative 
to the 2+ it does not exhibit the negative-?/ behavior 
(E f D (+) > {E f D (0),E f D (2+)} for all Fermi energies) ob- 
served in GqWq. In addition G W stabilizes the nega- 
tive charge state for the split<110> interstitial. If, on 
the other hand, the configuration with the lowest energy 
is considered, LDA and Go Wo superficially give a more 
similar picture: the tetrahedral 2+ state is stable for 60- 
70% of the respective band gaps. While LDA then gives 
preference to the neutral split<110> for larger Fermi lev- 
els, the Go Wo corrections marginally stabilize the neutral 
hex configuration, in agreement with the earlier DMC 
calculations 9 . The actual energies and transition levels 
between LDA and GoWo, however, differ appreciably. 

Every point at which two lines in Fig. [3] cross corre- 
sponds to a charge-state transition level e q / q >. Bracht 
et al. have recently determined these for the silicon 
self-interstitial in high temperature diffusion experiments 
[I]. They identified two levels, at ss 0.1-0.2 eV and at 
w 0.4 eV above the valence-band maximum, that they 
ascribed to e /+ an( l £ +/2+, respectively. These would 
most closely correspond to the GoWo-corrected charge- 
state transition levels £ / + =0.09 eV and e + / 2 +=0.58 eV 
for the hexagonal configuration or e /+=0-05 eV and 
e + /2+=0.50 eV for the split<110>, while those of the 
C^ v configuration are noticeably different (0.62 eV and 
1.24 eV). Although lowest in formation energy and there- 
fore highest in concentration, the C^ v 2+ configuration 
(which is identical to tet 2+) most likely plays a negligible 
role in the diffusion experiments, since its diffusion would 



have to proceed through a hexagonal site. The activation 
energy for this process (formation energy + energy bar- 
rier at the experimental situation of a Fermi level close 
to the middle of the band gap [I]) would thus be consid- 
erably larger than the activation energy for diffusion pro- 
cesses involving the other configurations. Refinements in 
the diffusion models (e.g. inclusion of multiple configura- 
tions and charge-state dependent diffusion barriers) may 
be able to clarify the role of the tet 2+ configuration in 
future experimental studies. 
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